Using modified Gaussian distribution to study the physical properties of one and 

two-component ultracold atoms 
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Gaussian distribution is commonly used as a good approximation to study the trapped one- 
component Bose-condensed atoms with relatively small nonlinear effect. It is not adequate in dealing 
with the one-component system of large nonlinear effect, nor the two-component system where 
phase separation exists. We propose a modified Gaussian distribution which is more effective when 
dealing with the one-component system with relatively large nonlinear terms as well as the two- 
component system. The modified Gaussian is also used to study the breathing modes of the two- 
component system, which shows a drastic change in the mode dispersion at the occurrence of the 
phase separation. The results obtained are in agreement with other numerical results. 

PACS numbers: 03.75.Hh, 32.80.Pj, 03.65.-w 
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I. INTRODUCTION 

Single-component Bosc-Einstein condensate (BEC) in 
trapped atoms has been observed more than ten years 
ago. Since then, in particular in recent years, laser cool- 
ing technique has been extended to trap atoms of multi- 
ple components. The basic technique is called "sympa- 
thetic cooling" that cools the second kind of atoms by 
means of the first kind pre-cooled atoms. Various com- 
binations of two-component atomic gas have been mixed 
and cooled. The first was done on the mixture of the 
same kind of bosonic atoms with different hyperfine spin 
states. For instance, the combinations of 87 Rb|2,2) and 
87 Rb|l,-l) states Q, and of 87 Rb|2, 1) and 87 Rb|l,-l) 
states [2|- Soon after, the experiment has been extended 
to cool the mixture of different kinds of bosons. This in- 
cludes the mixture of 87 Rb and 41 K atoms Q and 133 Cs 
and 87 Rb atoms [J] . The mixtures of bosons and fermions 
have also been achieved, such as the combination of 7 Li 
and 6 Li atoms [10, 6 Li and 23 Na atoms 0, and 40 K 
and 87 Rb atoms [8|. The fermion-fermion mixtures have 
also been accomplished, such as the combinations of two 
different spin states in 40 K atoms @ and in 6 Li atoms 
jTol . [lll | . The latter has led to the formation of bosonic 
molecules using Feshbach resonance and has advanced 
the observation of BEC-BCS crossover [H [H Ei . More 
recently, imbalanced spin population in fermion-fermion 
mixture has also been created and boosted researches on 
these systems [H, HU . 

In the literature, there exists theoretical works on equi- 
librium and dynamical properties of the two-component 
cold atoms, includin g th e investigations on the boson- 
boson mixtures fl7l. Tl8l. [l9l . I20I , l2l| and the fermion- 
fermion mixtures |22L 1 2 3^ For two-component boson- 
boson mixtures, in which both components are Bose 
condensed, there are two major theoretical approaches. 
One is concerned with the Thomas-Fermi approxima- 
tion (TEA), in which kinetic-energy terms are ignored. 
The other starts with the coupled Gross-Pitavcskii (GP) 
equations, which arc then solved numerically. In this pa- 



per, we attempt an alternative approach to investigate 
the two-component system based on variational method. 
Through the introduction of a modified Gaussain trial 
function, one is able to study the condition of the sys- 
tem being in the state of phase separation or misciblc. 
The breathing modes of the system is also studied and 
shown to be intimately connected to their ground-state 
behavior. As a matter of fact, the modified Gaussian dis- 
tribution is clearly shown to be a better description than 
the simple Gaussian even to the one-component trapped 
BEC, especially when the nonlinear effect is more im- 
portant. Variational method in adjunction with the in- 
troduced trial modified Gaussian function can be easily 
generalized to study the two-component boson-fcrmion 
and fermion-fermion mixtures. The latter includes the 
system of imbalanced spin population of current inter- 
est. 



The paper is organized as follows. In Sec. HH wc first 
give the Hamiltonian of a two-component boson-boson 
mixture. Coupled Gross-Pitaveskii energy functionals 
are derived for the system assuming that both compo- 
nents are fully Bose condensed. A subsection is devoted 
to sketch the variational method. In Sec. IIIII we intro- 
duce a trial "modified Gaussian function" whose ground- 
state behaviors to the one-component Bose system are 
compared in details with those of the simple Gaussian 
distribution and the TFA. In Sec. HVi the trial mod- 
ified Gaussian function are applied to study the equi- 
librium and dynamical properties of the two-component 
Bose system, and in Sec. [V] a brief conclusion is given. 
In the current context, we focus on the trapped system of 
boson-boson mixture. Similar approaches to the fcrmion- 
boson and fermion-fermion mixtures will be possible fu- 
ture works. 
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II. FORMALISM AND METHODOLOGY 



B. Variational Method 



A. Coupled GP Energy Functionals 

For a trapped two-component boson-boson mixture, 
the Hamiltonian can be given by 



H = H 1 +H 2 + H 



12- 



(1) 



where 



Hi = J dv4>\{v) 



and 



h 2 v 2 

2rrii 



+ Vi(r)-fXi 



drdr/4(r)4(v')5(r - r')^(r')^(r)(2) 



One major goal of this paper is to apply the variational 
method to investigate the dynamics of the coupled two- 
component ultracold atoms. The key to the variational 
approach is to first look for a trial wave function $i(r) 
for each component i, which is reasonable to the ground 
state of the system. With these trial wave functions, 
one then incorporates appropriate dynamical variables 
[$i(r) — > $i(r,t)], which are sensible to the dynamics 
under studies. The next step is to derive the correspond- 
ing Lagrangian: L = J drT — E using $j(r, t). Here 



T 



E 

t=l,2 



ih 
T 



<lv! 



at 



at 



(T) 



and 



H u = g l2 I / drdr'</»t(r)4(r')<5(r-r')</>2(r')^i(r!)3) 



Here i = 1,2 denote the species i of the mixture, </>j(r) 
and <^>j (r) are creation and annihilation operators on the 
particle of species i at r, and /ij denote the chemical po- 
tential of the species i. The entire system is trapped 
in an isotropic magnetic potential Vi(r) 



^mitofr 2 re- 



spectively for each component i. It is assumed that the 
interaction is dominated by the s-wave approximation in 
the dilute limit. Consequently the interaction coupling 
constants in (|2|) and ([3]) can be written in a unified way 



Airfra, 



9ij 



(4) 



E=(H}=E 1 +E 2 + V 12 , 



(8) 



where E±, E 2 , and Vi 2 are calculated based on (J5]) and 
([6]) with ^(r) being replaced by $i(r, t) now. Finally 
a set of coupled dynamical equations of motion can be 
derived through the Euler-Lagrange equations: 



dL 

dj3 



d dL 
di~df3' 



(9) 



where (3 is any one of the dynamical variables chosen. 
The corresponding dynamics, in particular the collective 
mode dispersion relation, can be obtained by solving the 
roots of the dynamical matrix (diagonalization) . 



where Ojj is the scattering length for a particle in species 
i with a particle in species j. The "mass" my = 
2m,mj / (rrij + rrij) with rm the bare mass of the particle 
in species i. For simplicity, we shall denote gn = g\, 
.922 = 52 and an = ai; 0,22 = »2 afterwards. 

At sufficiently low temperatures, consider that both 
components are fully Bose condensed, the thermal av- 
erage of Hi in @ then reduces to the Gross-Pitavcskii 
(GP) energy functional: 



(H- 



dr$*(r) 



h 2 X7 2 1 

■- + -m^Xr- -m 

Imi 2 
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fj< 



<f 4 (r)| 2 *i(r), 



(5) 



where $i = (<f>i). Similarly the thermal average of H12 
in j3]) reduces to 



Vv. 



(H 12 )=g 12 J < ir|* 1 (r)| 2 |* 2 (r) 



(0) 



III. TRIAL WAVE FUNCTION 

Before a proper trial wave function is introduced for 
each species in the trapped two-component cold atoms, 
it is useful to first examine the validity of the simple 
Gaussian distribution. 



A. Remarks: Gaussian Distribution 

When temperature approaches zero, single-component 
trapped Bose condensate (isotropic case) is governed by 
the GP energy functional given in Eq. ([5]). [In this case 
the index i in is redundant and is thus removed.] 
When the coupling g is infinitesimally small, the ground 
state of the system is perfectly a Gaussian distribution: 



$(r) = Cexp 



-Ar 2 ) 



(10) 



Gaussian function is "non-perturbative" and has the ad- 
vantage of easily integrated over the space. When g 
starts to increase, however, a Gaussian may no longer 
be a good one to describe the system. In the limit of 
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large g, Gaussian distribution fails and the most common 
approximation to the ground state of the system is the 
so-called "Thomas- Fermi approximation" (TFA), which 
ignores the kinetic term in ((5|) . The validity of the Gaus- 
sian distribution upon the increase of g can be sorted out 
as follows. First, let's rescale the GP energy functio nal: 
E/(Hlj) -> E, tuj -> t, r/i ->■ r with I = ^h/muj be- 
ing the length scale throughout this paper. Next, divide 
the energy functional by TV (total number of particles), 
Eq. (J51) for the current one-component system then be- 



E/N = J dr$* 



I r 2 + ^AT ff |$|2 

2 2i 



$. (11) 



By writing this way, it is seen clearly that in addition to 
g, both the magnetic trap potential (through £) and TV 
are also coupled to the nonlinear term. As a consequence, 
whether a Gaussian is a good approximation will depend 
on the magnitude of the front factor for the nonlinear 
term, mNg/2£ cx TV 'gvr?l 2 'y/uj. Provided that g is fixed, 
the smaller the value of Ny/oJ, the better the Gaussian 
distribution is for the same atoms. 

In current experiments performed on the single- 
component system, nonlinear effect is typically not small. 
Application of the Gaussian distribution to a real system 
thus always poses a question mark. Moreover, in a sys- 
tem of two-component mixture, the phenomenon "phase 
separation" could exist when the interspecies coupling 
<?i2 is strong. In this case, using a Gaussian to describe 
the distribution of the individual component will be fun- 
damentally ineffective. 



B. Modified Gaussian Distribution 

To better describe the system and for capturing the ad- 
vantage of easy integration of the Gaussian function, we 
propose the following "modified Gaussian distribution" : 



$(r) = cjexp 



A{r~r f 



cxp 



A(r + r ) 2 



C'exp ( — ) cosh(Ar r), 



(12) 



which is the superposition of two Gaussian functions, 
with centers deviated from the origin by ro and — ?*o re- 
spectively. As written explicitly in the second line of P^|) . 
the resulting function is a Gaussian function modified 
(multiplied) by a hypercosine function. In (TT2")) . C and 
C = 2Cexp(— Atq/2) are constants for normalization. 
A and ro are parameters corresponding to the amplitude 
and overall density profile of the atom cloud. For a given 
C, the values of A and ro are determined upon minimiz- 
ing the ground-state energy E of the system. In Fig. [TJ 
ro dependence of the normalized density profile |$(r)| 2 r 2 
is shown with a fixed A. As can be seen clearly, the 
maximum of |<f>(r)| 2 r 2 has a tendency to move outwards 




r„ = Q.4 
r„ = 0.8 
r„ = 1.2 
r„=1.6 
r = 2.0 
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FIG. 1: The normalized density profile, |<J?(r)| 2 r 2 , with $ 
given by the modified Gaussian function in (|12|l . Here A = 0.8 
is fixed for all cases. 



when ro is increased. When ro is large, the maximum 
could move significantly away from the origin. Having 
this feature, which is important when phase separation 
occurs, the modified Gaussian function is considered to 
be more effective for a two-component mixture. 

Compared to the simple Gaussian distribution 
[Eq. (fT0|) ]. the modified Gaussian consists one more pa- 
rameter (ro), that allows for a better description of the 
ground-state equilibrium properties of the system. As 
an illustration, in Fig. [2] we plot the ground-state energy 
E [Eq. (fTTj) ] as a function of the scattering length a for 
the modified Gaussian distribution. The parameters are 
taken to be: u = 2ir x 160Hz and TV = 10 4 . Similar re- 
sults for the simple Gaussian distribution and the TFA 
are also shown for comparison. In the small-a regime 
(small nonlinear effect), the curve of modified Gaussian 
merges to the one of simple Gaussian, as it must be. 
Gaussian distribution is exact in the a — > limit. The 
ground-state energy E is seen to be much underestimated 
for TFA at small a. In the intermediate a (= 2 ~ 5) 
regime, nonlinear effect becomes more and more impor- 
tant, the curve of modified Gaussian starts to deviate 
(becomes lower) from the one of a simple Gaussian. This 
indicates that the modified Gaussian distribution is more 
favored by the system so long as the ground-state energy 
is concerned. It is worth noting that most current exper- 
iments lie in this intermediate a regime. In the regime of 
large a, the results of TFA should be the most accurate 
ones. 

In Fig. [31 spatial distributions (density profile) of the 
three cases (Gaussian, modified Gaussian, and TFA) are 
compared for different total numbers of particle TV. The 
case of larger TV corresponds to the case of larger non- 
linear effect. In frame (a) of smaller TV, the curve of 
modified Gaussian is closer to the one of simple Gaus- 
sian - with the maximum at the origin. When N is larger 
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FIG. 2: Ground-state energy E divided by hio (u = trap 
frequency) versus the scattering length a for single-component 
condensed bosons. The curves are for Gaussian (G), modified 
Gaussian (MG), and TFA. The parameters are lj — 2n x 
160Hz and N = 10 4 . 



[frame (b) & (c)] , the curves of modified Gaussian behave 
closer to the ones of TFA, although the maxima may de- 
viate from the origin. As far as a real distribution is 
concerned, this feature is somewhat defective. It should 
be emphasized, however, that modified Gaussian has an 
overall distribution which is closer to the one of TFA 
and is more accurate for large nonlinear effect. More- 
over, modified Gaussian inherits the advantage of easy 
integration of simple Gaussian, which makes variational 
method feasible for the problems. 

In next section, equilibrium property and dynamics of 
a two-component system arc studied. Breathing modes of 
such system will be shown to depend crucially on whether 
the system is miscible or phase-separated. 



IV. TWO-COMPONENT MIXTURE 

The analyses in previous section are limited to single- 
component systems. In the case of two-component mix- 
ture, the system could undergo a phase transition from a 
miscible state to a phase-separated state when the inter- 
species coupling is strong. For current inhomogencous 
but isotropic system, when phase separation exists, one 
component of the atoms will have maximum density de- 
viated from the origin, while the density of the other 
component remains maximum at the origin. As pointed 
out before, a simple Gaussian is ineffective to describe 
such phenomenon. The present section is devoted to 
study how the phase separation occurs when the modified 
Gaussian distributions are in effect. At the same time, 
the dynamics (breathing modes) of the two-component 
system will also be studied starting from the same mod- 
ified Gaussian distribution. 
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FIG. 3: The spatial dependence of the density profile |$(r)| 2 
for a one-component system. Three kinds of distribution 
(Gaussian, modified Gaussian, and TFA) are compared. 
Frame (a)-(c) are for number of atoms N = 200, 10 3 , and 
10 4 respectively. 



A. Equilibrium Property 

The ground-state properties of the two-component 
bosonic atoms are first studied. Considering that both 
components are described by the modified Gaussian dis- 
tribution 



$i(r) = C;exp (-^A,r 2 



cosh^r^r), (13) 



where i = 1, 2 and the normalization factor Cj is depen- 
dent of the total number N t of individual component i. 
The ground-state energy E of the system can be calcu- 
lated by substitution of flT3J) into Eqs. (HJ)-® and ©. 
Given three coupling strength oi, 02, and 012, and the 
particle number N\ and N2, final result of E can be ob- 
tained by minimizing E against the variation of A\ and 
A2, and r± and r2- In Fig. |H the results of E are shown 
as a function of the inter-component scattering length 
ai 2 . Two cases of equal population N\ = N 2 = N = 10 3 
and N = 10 4 are presented. The slopes, dEjdayi vs. 
a 12 are also shown. For easy comparison, throughout 
this section theparameters are taken to be the same as 
those in Ref. [nj|. Based on the experiment of Na-Rb 
mixture, the scattering lengths of Na and Rb atoms are 
chosen to be 3 nm and 6 nm respectively. Regarding the 
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FIG. 4: The dependence of inter-component scattering length 
ai2 on the ground-state energy E, divided by Nhuai, of a two- 
component Boson system. Modified Gaussian distribution is 
undertaken. The slopes of dE/dai2 vs. ai2 are also shown 
[E' = E/(NhuJi)]. Left column corresponds to Ni — N2 = 
N = 10 3 , while right column corresponds to TV = 10 4 . See 
text for the values of other parameters. 



external magnetic trap, the corresponding trap frequen- 
cies are lo\ = 2tt x 310Hz for Na and u>2 = 27r x lQOHz 
for Rb respectively. In both plots in Fig. 0] (especially 
dE/dai2 vs. ai 2 ), one sees clear evidence that the sys- 
tem undergoes a transition from one phase to another at 
a 12 ~ 2 — 2.5 when N = 10 4 and at an ^4—5 when 
N = 10 3 . 

Fig. [5] shows how this transition behaves in the den- 
sity profiles of the two-component system. The figures 
are plotted in a variety of interspecies scattering lengths 
(012) with fixed a% = 3 nm and ai = 6 nm. In the case of 
N = 10 4 , it is found that both |$i(r)| 2 and |$ 2 (r)| 2 nave 
maximum away from the center when ayi = 0. It has to 
be emphasized again, that modified Gaussian distribu- 
tion has a feature that the maximum density could de- 
viate from the origin due to large nonlinear effect. Even 
though ai2 = 0, the individual nonlinear effect is still 
large, which lead to such "unphysical" result. As men- 
tioned before, it is the fact that modified Gaussian has an 
overall distribution more close to TFA, that one has to 
look into. When becomes larger (012 ~ 2 — 2.5), the 
maximum of |$i(r)| 2 is pushed significantly away from 
the center, while the maximum of | $2(7) 1 2 is right at the 
center - the pheomenon of phase separation. 

To see this phenomenon more clearly, in Fig. [5] we plot 
the value of parameter t\ (for the modified Gaussian) 
versus a\i. It is reminded that r\ is referred to the Na 
component in Fig. 5. In the case of N = 10 4 , there 
is a sudden change for n at a\2 ~ 2 — 2.5. This is in 
great consistence with the ground-state energy shown in 
Fig. [4] where a sharp change also occurs at this range 
of a\2- In view of the original definition of r\ in the 
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FIG. 5: The individual density profile plotted in a variety 
of the interspecies scattering length 012 for a two-component 
system. From the top to the bottom are for N\ = JV2 = N = 
10 4 , 10 3 and 200 respectively. Note that in x axis, r is scaled 

to £2 — ^h/m,2tjJ2- 



first line of Eq. (fT!?)) (i.e., ro there), a sudden change of 
r\ means a sudden change of the coordinate of density 
maximum (see also Fig. 1). This sudden change indeed 
corresponds to the transition of "phase separation" . In 
the case of N = 10 3 , we also see a consistent picture 
among the results in Figs.HHS] A phase separation occurs 
at a\2 ~ 4 — 5. When = 200, the phenomenon is less 
obvious. However, one can still see a phase separation 
that appears at a\i ~ 9 — 10. 



B. Breathing Modes 

In this section, we investigate how the phase separa- 
tion affects the dynamics of the two-component system. 
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FIG. 6: The value of the parameter ri divided by £2 (r'i is 
referred to Na component in Fig. 5) for the modified Gaussian 
plotted as a function of the interspecies scattering length a\2- 

It was shown in Ref. that for a two-component mix- 
ture, the collective mode dispersion will have a drastic 
change near the occurrence of phase separation. Since 
the property of collective mode is intimately connected to 
the equilibrium property of the system, the sharp transi- 
tion in the ground-state energy of the system could imply 
a drastic collective mode dispersion change in the same 
regime. 

In the current context, we will focus on the breathing 
modes of the system. One first extends the modified 
Gaussian distribution for each com pon ent [Eq. t| 13|>] to 
include the dynamical variables [24l, |25|: 

§i(r,*) = C l exp^~^A l [l + e l (t) + ie' l {t)}r 2 ^ 

x cosh^r-ir). (14) 

Here £j and e\ correspond to the fluctuations of local 
amplitude and local phase of the atom cloud along the 
radial r direction. Ai and are determined earlier 
through minimizing the ground-state energy of the sys- 
tem (see previous subsection). Since we are interested 
in the breathing modes of the lowest energy and of only 
one radial node (n = 1), and are coupled only to 
the r 2 term in the Gaussain function in Q14p . One may 
introduce extra dynamical variables to the hypcrcosinc 
function in (fl4|) . that could lead to more branches of ra- 
dial modes. Some of them will have more than one node 
(n > 2). 

Substitution of the above dynamical wave function into 
Eqs. Q and ©, one can derive the coupled equations of 
motion [using ©] for Si and e' i: and then solve a corre- 
sponding dynamical matrix to obtain the dispersions of 
relations for the breathing modes. In Fig. [3 we present 
the breathing mode dispersion as a function of the intcr- 




a |2 (nm) 

FIG. 7: (Color online) Based on the variational method (V), 
in-phase and out-of-phase breathing mode dispersions are 
plotted as a function of ai2 for a two-component boson sys- 
tem. In each frame, first turning point in lower curves corre- 
sponds to where the system starts to be phase separated. The 
second turning point corresponds to the exchange of the in- 
phase and out-of-phase modes. In the left frame of N = 10 4 , 
the results given by numerical calculation [19] (N) are com- 
pared. The two crosses at a\2 — correspond to the breathing 
modes of individual component 1 and 2 (decoupled regime). 



component scattering length an. Again, all parameters 
are the same as before. In the case of N = 10 4 , the results 
given by numerical calculation in Ref. (l9l | is included for 
comparison. 

Much information is revealed in Fig. [7] Taking the 
N = 10 4 case as an example, when interspecies scattering 
length a\2 is negative, the system is miscible. In-phase 
breathing mode has lower energy and is more easily ex- 
cited. When ai2 is positive, the results are divided into 
three regimes. When < a\% < 2.5, the system is mis- 
cible and thus out-of-phase breathing mode has lower 
energy and more easily gets excited. When the system 
starts to phase separate but most parts of the two com- 
ponents are still overlapping (2.5 < a.12 4), the energy 
of the out-of-phase breathing mode remains lower than 
that of in-phase mode . When the two components are 
well separated (not much overlap) at a.12 > 4 (see also 
Fig. 5), in-phase breathing mode turns out to have lower 
energy again. The switch of the in-phase and out-of- 
phase breathing modes in the above latter two cases can 
be understood alternatively as follows. When phase sepa- 
ration occurs but the two components are still very much 
overlapped, out-of-phase mode is more easily excited be- 
cause it corresponds to the decrease of the overlap (and 
the energy). However, when the two components are well 
separated, in-phase mode is in turn more easily excited 
because in this case out-of-phase mode simply increases 
the overlap (and the energy). 

As seen in Fig. [TJa), there are two turning points in 
the lower dispersion curve. The first turning point is the 
signature of the system being phase separated. While 
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the second turning point corresponds to the exchange 
of the in-phase and out-of-phase breathing modes in the 
regime where the two components are well separated. In 
comparison with the numerical results given by Ref. [19] 
(where only the lowest-energy mode is given), our result 
of lowest-energy mode matches theirs quantitatively be- 
fore phase separation arises. After the phase separation 
occurs, our result is qualitatively similar to their. (There 
is no second turning point seen in Ref. [19] for the range 
plotted.) It is worth noting that the two points of a 12 = 
correspond actually to the breathing modes of two decou- 
pled components (lower-frequency one is for Rb, while 
higher- frequency one is for Na). From the goodness of 
the match of our results to those of Ref. [19], it strongly 
suggests that modify Gaussian is also a good function to 
describe the dynamics of a one-component system when 
the nonlinear effect is important. 

V. CONCLUSIONS 

In conclusion, the density profile, the ground-state 
energy, and the breathing modes of a two-component 



Bose condensed system are studied using a variational 
method. We propose a modified Gaussian distribution 
function which is shown to be more effective for a two- 
component system. While the modified Gaussian func- 
tion reserves the advantage of easy integration with the 
simple Gaussian, it actually gives a better description 
even to the equilibrium and dynamical properties of a 
single-component system. Our approach can be gener- 
alized to study the vortices of a two-component system, 
and possibly the system of unbalanced spin population. 
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